In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
from sklearn import datasets, linear_model
from scipy.stats import linregress
from numpy.linalg import inv
import statsmodels.discrete.discrete_model as sm
from numpy.linalg import eigvals
import numpy.linalg as LA
%matplotlib inline
def set_data(p, x):
temp = x.flatten()
n = len(temp[p:])
x_T = temp[p:].reshape((n, 1))
X_p = np.ones((n, p + 1))
for i in range(1, p + 1):
X_p[:, i] = temp[i - 1: i - 1 + n]
return X_p, x_T
In [2]:
data = sio.loadmat('Tut3_file1.mat')
DLPFC = data['DLPFC']
x = DLPFC[0, :]
plt.plot(DLPFC.T)
Out[2]:
In [43]:
x = DLPFC[0]
X_p, x_T = set_data(3, x)
In [44]:
res = LA.lstsq(X_p, x_T)
In [69]:
eps = x_T - np.dot(X_p, res[0])
XpXp = np.diag(LA.inv(np.dot(X_p.T, X_p)))
sdeps = np.std(eps, axis=0, ddof=0)
res[0].flatten() / (sdeps * np.sqrt(XpXp))
err = (sdeps * np.sqrt(XpXp))
Out[69]:
The t-statistic computes how distant our parameters are from the zero hypothesis in terms of standard errors. So for $a_0$ and $a_1$ we cannot refute the zero hypothesis, but in case of $a_2$ and $a_3$ it is highly unlikely that our that comes from random fluctuations.
In [53]:
import statsmodels.api as sm
ols = sm.OLS(x_T, X_p).fit()
ols.summary()
Out[53]:
The likelihood function is: $ f (Y) = (\frac{1}{\sqrt{2 \pi} \sigma} \exp{- \frac{1}{2 \sigma^2}\sum^n_{i=1}(Y_i - a_0 X_{i1} - ... - a_4)}) $, so the log likelihood will be: $-n\ln(\sqrt{2 \pi} \sigma) - \frac{1}{2 \sigma^2}\sum^n_{i=1}(Y_i - a_0 X_{i1} - ... - a_4)$. The result in our case can be seen in the summary of the previous cell.
In [70]:
ols.params, res[0].flatten()
ols.HC1_se,
Out[70]:
In [18]:
X = np.ones((5, len(x)))
X[1:3, :] = DLPFC
X[3:, :] = data['Parietal']
X_T = X[1:, 1:].T.reshape((359, 4))
A = np.zeros((4, 5))
tv = np.zeros((4, 5))
for i in range(1, 5):
X_T = X[i, 1:].T.reshape((359, 1))
X_p = X[:, :-1].T
ols = sm.OLS(X_T, X_p).fit()
A[i-1, :] = ols.params
tv[i-1, :] = ols.pvalues
We can see on the off diagonal terms are much smaller than the terms on diagonal. That shows that there is not as much influence between features as on the feature itself.
In [20]:
A[:, :]
Out[20]:
In [29]:
X_T = X[1:, 1:].T.reshape((359, 4))
res = LA.lstsq(X_p, X_T)
For checking if it si astationaty process we will see if the eignvalues are in the unitcircle. They are.
In [33]:
res[0].shape
Out[33]:
The loglikelihood of this model will be the
In [36]:
eps = X_T - np.dot(X_p, res[0])
sdeps = np.std(eps, axis=0)
In [39]:
Out[39]:
In [41]:
XpXp = LA.inv(np.dot(X_p.T, X_p))
In [ ]: